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ABSTRACT 

There is compelling evidence that most -if not all- galaxies harbour a super-massive 
black hole (SMBH) at their nucleus, hence binaries of these massive objects are 
an inevitable product of the hierarchical evolution of structures in the universe, 
and represent an important but thus-far elusive phase of galaxy evolution. Gas 
accretion via a circumbinary disc is thought to be important for the dynamical 
evolution of SMBH binaries, as well as in producing luminous emission that can be 
used to infer their properties. One plausible source of the gaseous fuel is clumps 
of gas formed due to turbulence and gravitational instabilities in the interstellar 
medium, that later fall toward and interact with the binary. In this context, we model 
numerically the evolution of turbulent clouds in near-radial infall onto equal-mass 
SMBH binaries, using a modified version of the SPH code gadget-3. We present a 
total of 12 simulations that explore different possible pericentre distances and relative 
inclinations, and show that the formation of circumbinary discs and discs around each 
SMBH (‘mini-discs’) depend on those parameters. We also study the dynamics of the 
formed discs, and the variability of the feeding rate onto the SMBHs in the different 
configurations. 

Key words: accretion, accretion discs, black hole physics, hydrodynamics, galaxies: 
evolution, galaxies: nuclei 


1 INTRODUCTION 

The idea that super-massive black holes (SMBHs) reside 
in the nuclei of at least the more massive galaxies is 
currently well established (Richstone et al. 1998; Kormendy 
& Ho 2013). Additionally, according to our current structure 
formation paradigm, galaxies often interact and merge 
following the hierarchical growth of their parent dark matter 
halos. Then, in the aftermath of a major galaxy merger, 
we expect that the merger product will contain a pair of 
SMBHs. Dynamical interaction with stars and gas drive 
the SMBHs towards the centre of the new galaxy on short 
time-scales, where they will form a gravitationally bound 
black hole binary (BHB; Begelman et al. 1980; Milosavljevic 
& Merritt 2001). The further evolution of the BHB will 
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depend on the content and distribution of both gas and 
stars in the galactic nucleus (see e.g. Sesana 2010; Khan 
et al. 2013; Holley-Bockelmann & Khan 2015). 

In gas-rich systems, the BHB evolution is likely to be 
driven primarily by gas. Many theoretical and numerical 
studies have focused on the orbital evolution of a sub-parsec 
binary surrounded by a gaseous circumbinary disc (Ivanov 
et al. 1999; Gould & Rix 2000; Armitage & Natarajan 
2005; Cuadra et al. 2009; Haiman et al. 2009; Lodato 
et al. 2009; Nixon et al. 2011a; Roedig et al. 2011; Kocsis 
et al. 2012; Nixon 2012; Roedig et al. 2012; Amaro-Seoane 
et al. 2013b; Nixon et al. 2013; Roedig & Sesana 2014) 
and the electromagnetic signatures that such a system 
would produce (Artymowicz & Lubow 1996; Milosavljevic & 
Phinney 2005; Tanaka & Menou 2010; Tanaka et al. 2010, 
2012; Sesana et al. 2012; D’Orazio et al. 2013; Tanaka 2013; 
Tanaka & Haiman 2013; Roedig et al. 2014; Brem et al. 
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2014; Farris et al. 2015). However, the exact mechanism that 
would form such discs is still unclear, as it critically depends 
on the fuelling processes onto those scales. 

Observational (e.g. Sanders & Mirabel 1996) and 
numerical studies (e.g. Mayer et al. 2007) show an 
abundance of dense gas in the nuclei of merged galaxies. 
However, the wide range of physical scales involved makes 
the fuelling a very complex process - the gas needs to 
lose its angular momentum efficiently to be transported 
from galactic scales down to the nuclear region (Hopkins 
& Quataert 2011). An important factor controlling the 
fuelling is the density distribution of the gas. Turbulence 
and gravitational instabilities in the interstellar medium 
can trigger the formation of gas clumps, which will be the 
seeds for molecular clouds (Agertz et al. 2009). Relaxation 
processes might then randomise the orbits of the clouds, 
producing discrete events of “ballistic” accretion onto a 
central SMBH, almost unaffected by hydrodynamical drag. 

The large spread in angular momentum achieved with 
turbulence is likely to -at least partially- randomise the 
direction of the accretion events even if the bulge has 
a net rotation (Hobbs et al. 2011). A scenario in which 
SMBHs evolve by accreting individual gas clouds falling 
from uncorrelated directions has been proposed by King & 
Pringle (2006). In later publications in the context of BHBs, 
Nixon and collaborators have shown that this scenario can 
lead to the formation of counter-rotating discs, facilitating 
the binary ffnal coalescence (Nixon et al. 2011a,b; Dunhill 
et al. 2014; Aly et al. 2015). In a somewhat different 
application, the feeding of SMBHs through individual cloud 
infall events with different degree of angular momenta 
randomisation has been proposed as a viable scenario to 
reproduce current measurement of SMBH spins (Dotti et al. 
2013; Sesana et al. 2014). 

A perhaps more concrete example of this kind of 
accretion events is the putative cloud that resulted in the 
unusual distribution of stars orbiting our Galaxy’s SMBH. 
Several numerical studies have shown that portions of a 
near-radial gas cloud infall can be captured by a SMBH to 
form one or more eccentric discs that eventually fragment 
to form stars (Bonnell & Rice 2008; Hobbs & Nayakshin 
2009; Alig et al. 2011; Mapelli et al. 2012; Lucas et al. 
2013), roughly reproducing the observed stellar distribution 
(Paumard et al. 2006; Lu et al. 2009). 

Working on the hypothesis that infalling clouds are 
common in post-merger galactic nuclei hosting BHBs, 
we study the formation and early evolution of discs in 
these events. Recently, Dunhill et al. (2014) showed that 
misaligned clouds infalling with a large impact parameter 
onto BHBs result in prograde or retrograde circumbinary 
discs. In this study, we focus on cloud orbits with lower 
speciffc angular momentum, as well as higher inclinations 
relative to the BHB orbit, in order to investigate a wider 
variety of disc conffgurations and dynamics. 

The paper is organised as follows. In Section 2, we 
describe the details of our numerical model. Section 3 
discusses the output of the simulations and the formation of 
discs depending on the inclination and impact parameter. 
We detail in Section 4 the accretion on to the binary and 
each SMBH, in particular its variability. In Sections 5 and 6, 
we study the dynamics of the mini-discs and circumbinary 
discs, respectively. Finally, we discuss the implications of 


our results on the observability of these systems and future 
multi-messenger studies in Section 7, and summarise our 
study in Section 8. This paper is the first on a series where we 
will be reporting the results of numerical models of infalling 
clouds onto BHBs; in upcoming publications we will study 
the orbital evolution of the binary, detailed observational 
signatures from the formed discs, and the impact of magnetic 
fields, among other topics. 


2 THE NUMERICAL MODEL 

We model the interaction between the binary and 
cloud using a modified version of the smoothed particle 
hydrodynamics (SPH) code gadget-3, which is an improved 
version of the public code gadget-2 (Springel 2005). 

The binary is modelled by two equal-mass “sink” 
particles (see below) with initially circular, Keplerian orbits. 
The cloud is represented using approximately 4.2 x 10® 
equal-mass SPH particles, with a total mass 100 times 
smaller than the binary. 

The code units are such that the initial semi-major 
axis, mass and orbital period of the binary are unity. As 
a reference, this corresponds for one of our fiducial models 
presented in Section 7 to a lO^M© binary, with a semi-major 
axis of 0.2 pc and a period of ~ 8400 yrs. 

2.1 Initial conditions 

The cloud is initially spherical, with uniform density and 
a turbulent internal velocity held. The main effect of the 
turbulence is to provide support against the initial collapse 
due to the cloud’s self-gravity. The velocities are drawn 
from a Gaussian random distribution with power spectrum 
Pv{k) = (|Ffc|^) oc where k is the wave number of 

the velocity perturbation, and the power index is chosen to 
match the observed velocity dispersion profile of molecular 
clouds (Larson 1981). 

To set the internal velocity held we treat the cloud 
as incompressible (V • u = 0), which implies that we can 
represent the velocity held with a “vector potential" A. We 
sampled the components of A in the Fourier space using 
an equi-spaced lattice with 256® coordinates. To assign a 
velocity vector to each particle, we inverse Fourier transform 
A and we interpolate the obtain values between grid points. 
Finally, the velocity held is normalised such that the kinetic 
energy is equal to the absolute value of the potential energy, 
resulting in a cloud that is marginally unbound. 

The physical setup of the simulations is shown in 
Figure 1. It is based on the work by Bonnell & Rice (2008), 
but replacing the central single black hole with a binary. We 
place the cloud at a distance of 15 from the centre of mass of 
the binary, and we give it an initial velocity Vini such that it 
has an eccentric, bound orbit. The modulus of the initial 
velocity vector is constant in all our simulations (uini = 
0.5ubin, where Ubin = 0.5y^GM/a is the tangential velocity 
of each SMBH), but we change its direction, reflected in 
the angle ^vei- We model clouds with three different impact 
parameters; we choose ^vei = 0.197, 0.298 and 0.403 radians, 
so that the pericentre distances are 0.7, 1.5 and 3Rbin, 
respectively, where Rbin is the binary radius (i.e. half of the 
binary separation, a/2 = 0.5). 
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Figure 1. Initial setup of the simulations. The circle on the 
right represents the cloud, while the small black, solid circles 
are the SMBHs. In our set of 12 simulations we sample different 
inclinations of the cloud initial velocity (i.e. different ^vel) cind 
orientations of the binary orbit. 


Finally, since turbulence and relaxation might cause 
significant randomisation of the angular momenta of the 
gas clouds on parsec scales, we model clouds approaching 
with different relative orientations with respect to the binary 
orbital 

• “Aligned": the cloud starts in the plane of the BHB, 
with Vini lying on the same plane such that the cloud and 
the binary are co-rotating. 

• “Counter-aligned": same as aligned, but 

counter-rotating. 

• “Perpendicular edge-on": the cloud is initially in the 
same plane as the BHB, but the tangential component of 
Vini (r’ini sin^vei) IS perpendicular to that plane. 

• “Perpendicular face-on": the cloud starts in the plane 
perpendicular to that of the BHB, but the tangential 
component of Vini is parallel to that plane. 


2.2 Accretion 

The accretion recipe included in the standard version of 
gadget-3 uses the Bondi model to estimate the amount 
of mass that should be added to each SMBH and then 
determines, probabilistically, the corresponding accreted 
particles (Springel et al. 2005). This prescription for 
the black hole growth is usually applied in cosmological 
simulations (e.g. Planelles et al. 2014) because the scales 
where the SMBHs are accreting are well below the resolution 
limit. In contrast, in our numerical models we do resolve 
the Bondi-Hoyle-Littleton radius, thus we shall use a more 
deterministic recipe for accretion. 

In our simulations each black hole is represented by 
a “deterministic” sink particle. That means, it accretes all 
SPH particles satisfying some given conditions within a 
certain radius (Cuadra et al. 2006). In addition, this type 
of particle interacts with others through gravity, but not 
with hydrodynamical forces. In our model the SMBHs have 
a fixed accretion radius of rgink = 0 . 1 ; each particle crossing 
^sink is accreted if its kinetic energy is less than a fraction a 


^ Animations of our 12 simulations can be seen in http://www. 
astro.puc.cl/~fgarrido/animations 


of its potential energy^ (Dotti et al. 2009): 

^kin < Q:|^pot|. (1) 

Since our accretion radius is very large compared to the 
Schwarzschild radius^ (Rsch), this condition is necessary 
to avoid non-bound particles being added to the SMBHs. 
We adopt a = 1 throughout all the simulations presented 
here, meaning that all bound particles within the accretion 
radius are added to the corresponding SMBH. We did 
test a = 0.5 and also rgink = 0.05, but the results were 
virtually unchanged, thus we kept the original values to 
save computational time. We discuss how the accretion 
rate measured onto the sink particles can be interpreted in 
Section 4. 

The properties of the accreted particles, such as time, 
ID, mass, position and velocity, are stored once they 
are added to the SMBH (Dotti et al. 2010). Using this 
information we compute the total angular momentum of the 
gas accreted between outputs in the reference frame of the 
corresponding SMBH to explore the possibility of unresolved 
discs inside rgink in Section 5. 

2.3 Thermodynamics 

As we do not implement radiative cooling explicitly in our 
model, we use a barotropic equation of state, i.e. P = 
P{p). The functional form of the pressure is chosen to 
mimic the thermodynamics - in particular the temperature 
dependance with density - of star-forming gas. At low 
densities the cloud is initially optically thin to the thermal 
emission from dust grains, and the compressional heating 
rate by the collapse is much smaller than the cooling rate 
by the thermal radiation. The situation reverses at high 
densities, when the compressional heating overwhelms the 
radiative cooling rate and the gas is heated as collapses 
(Masunaga et al. 1998; Masunaga & Inutsuka 2000). 

The equation of state has the following form: 

P{p) = A{s)p'^, (2) 

where P is the pressure, p is the density, A(s) is the entropic 
function (Springel & Hernquist 2002 ) and 7 is the polytropic 
index that depends on the density as follows: 

7 = 1.0 for p ^ Pc, 

7 = 1.4 for p > pc, 

with Pc — 1.096 in code units. Note that the introduction 
of this two regime equation of state breaks the scale-free 
nature of our simulations. In Section 7.1 we discuss the 
interpretation of scaling this number to different physical 
units. The index 7 = 7/5 corresponds to an adiabatic regime 
for diatomic gas and it is the value found by Masunaga 
& Inutsuka (2000) that represent the thermal evolution of 
a collapsing molecular cloud in the density regime we are 
modelling. 

A consequence of using this equation of state is that 
we stop the collapse of the densest gas, avoiding excessively 

^ Both energies are computed in the reference frame of each 
SMBH. 

^ For a 10 ®Mq binary separated O.I pc, the sink radius would be 
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Figure 2. Column density map of an early snapshot after the 
initial conditions of one of the simulations. At this stage the 
cloud looks roughly the same for all orbital configurations, as 
its evolution is dominated by the initial turbulent velocity field. 
The SMBHs are represented with two black circles on the left. 

small time-steps that can stall the simulations. This type 
of equation of state is frequently used on hydrodynamical 
simulations of the evolution of turbulent clouds (e.g. Bate 
et al. 1995). In the case of our model, this simple treatment 
of the thermodynamics allows us to capture the global 
behaviour of the gas during the interaction with the 
binary without an explicit implementation of cooling and/or 
radiative transfer. As the overall gas dynamics during 
the early phases of the interaction is dominated by the 
gravitational potential of the binary, the results presented 
here will depend only weakly on the thermodynamics 
adopted. However, the long-term evolution of the gas will 
be likely dependent on the thermodynamics, e.g. the cooling 
rate will determine whether the gas fragments or not. 
However the modelling of those processes is beyond the 
scope of this paper. 


3 FORMATION OF DISCS 

In Fig. 2 we show an early snapshot of one of our simulations 
as a column density map. Due to the large initial distance 
between the cloud and the BHB, the evolution of the gas 
is initially dominated by the turbulent velocity field, which 
produces filaments in the cloud. In this snapshot we can 
already notice how the gas is stretched by the gravitational 
pull of the binary. At this stage, the gas evolves almost 
independently of the particular orbital configuration of 
the system, because the effects that the BHB quadrupole 
potential can have on the hydrodynamics of the cloud are 
negligible, and the differences between impact parameters 
are small. 

The study of the secular evolution of our systems, after 
the gas dynamics reaches a quasi-steady state, requires a 
considerable computing time that is not affordable with our 
standard configuration. We hence stop the simulations either 
once the transient effects of the cloud infall have stopped, 
or when the simulation stalls due to clump formation. We 
do however explore long-term effects with lower-resolution 
simulations in Section 6. 

We show the final state of our aligned simulations 
(models A) in the top row of Fig. 3. The first interaction 
of the gas with the binary is characterised by an efficient 


slingshot which pushes the remaining gas away from the 
system, shaping a tail. For the smallest impact parameter 
(top row, left panel of Fig. 3), the strong interaction only 
allows the formation of the so-called “mini-discs” around 
each SMBH, while increasing the impact parameter we 
allow more material to avoid the slingshot and settle in 
bound orbits around the binary, forming a circumbinary 
disc. The common feature of every impact parameter is the 
formation of mini-discs, appearing in the figure as black 
rings around either SMBH because of their high column 
density. The discs tend to be slightly misaligned with respect 
to the binary orbit. The analysis of the dynamics of these 
mini-discs is shown in Section 5. Only with this orientation 
we see the formation of extended, prominent and persistent 
mini-dises. For the other BHB-cloud orientations, the 
material captured by individual SMBHs has little angular 
momentum, thus falling directly within rgink- Nevertheless, 
in the higher resolution tests shown in Section 2.2 we do 
observe the formation of small and intermittent mini-discs, 
which indicates that mini-discs might form on smaller scales, 
but cannot be resolved with our standard resolution. We 
explore the possibility of unresolved mini-discs inside the 
sink radius in Section 5. 

The final state of our counter-aligned simulations 
(models CA) is shown in the upper middle row of Fig. 3. The 
interaction is very different because the typical gas velocity 
has the opposite direction relative to the orbital motion 
of the binary, cancelling most of the gas initial angular 
momentum. This enhances the accretion on to the BHB, 
as discussed in Section 4. The gas that remains bound after 
the interaction forms a very eccentric tail. For the three 
cases we observe the formation of a nearly circular, very 
narrow, counter-rotating circumbinary ring, arising as the 
material from the tail reaches the binary radius. The inner 
edge of these rings has a radius a from the center-of-mass 
(CoM) of the binary, which is expected due to the absence 
of resonances in a counter-rotating case (e.g. Nixon et al. 
201 la). The main qualitative difference observed on these 
three counter-aligned cases is the amount of gas clumps 
formed, as larger impact parameters result on more clumps. 
Most of these clumps form due to the compression of gas 
during the pericentre passage that allows them to become 
self-gravitating, and eventually to form stars. We discuss the 
observational implications of star formation in Section 7. 

We show the final state of our two perpendicular 
simulations in the lower rows of Figs. 3 (model PE, 
perpendicular edge-on and models PF, perpendicular 
face-on). As the typical gas velocity is perpendicular to that 
of the black holes, the slingshot is not very efficient, meaning 
that most of the gas is not pushed away from the binary. As 
we increase the impact parameter, we allow more material to 
settle on stable orbits. This gas forms a circumbinary disc, 
although completely misaligned with respect to the binary, 
keeping its initial angular momentum direction. In these 
configurations we also have regions of over-density which 
might lead to stellar formation, although they are located 
outside the region shown in the figures. 

In Table I we compile different relevant values we 
measure for the discs formed in our simulations. As noticed 
above, only the aligned cases show the formation of extended 
mini-discs, detectable at the resolution of our sink radius. 
These discs are very prominent and stable, and hence 
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Figure 3. Column density maps of our simulations. The different rows are the 4 inclinations modelled, while increasing the impact 
parameter from left to right, as indicated by the model name. The position of the SMBHs is indicated by the black crosses. Top 
row - aligned (A) simulations: The binary moves on the x-y plane, counter-clockwise. For every impact parameter dense “mini-discs" 
form around each SMBH. As the impact parameter increases, we observe the formation of circumbinary discs. Upper middle row - 
eounter-aligned (CA) simulations: The binary moves on the x-y plane, clockwise. There is no formation of mini-discs. For each impact 
parameter we observe the formation of an almost circular and retrograde circumbinary ring. On the right panel we see that the gas 
fragments due to the compression during the first passage and forms a large amount of clumps. Lower middle row - perpendieular 
edge-on (PE) simulations: The binary moves on the x-z plane. There is no formation of mini-discs. For the larger impact parameter 
the gas settles in bound orbits to form a circumbinary disc, which is clearly seen on the right panel. This disc is almost completely 
perpendicular to the orbit of the binary. Bottom row - perpendieular faee-on (PF) simulations: The binary moves on the y-z plane. 
Similar to the PE models, for the larger impact parameter the gas is able to form a perpendicular circumbinary disc. 
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Table 1. Properties of the discs formed in our simulations. In the model names the letter indicates the orbit orientation (A: aligned, 
CA: counter-aligned, PE: perpendicular edge-on, PF: perpendicular face-on) and the number the pericentre distance in units of binary 
radius, is the time at which we stop each simulation. Mi and r^^out (^ = 1,2) are the mini-disc masses and outer radii, respectively. 
Available M corresponds to the mass available to form a circumbinary disc, while < e > is the median eccentricity of that gas. 


Model 

^fin 

Resolved 

Ml 

M2 

'ri,out 

r2,out 

Noticeable 

Available M 

< e > 


(Pbin) 

mini-discs? 

(Mbin) 

(Mbin) 

(a) 

(a) 

Circumbinary? 

(Mbin) 


AO.7 

22.9 

YES 

2.6 X 10-4 

1.4 X 10-4 

0.34 

0.44 

NO 

2.9 X 10-® 

0.89 

A1.5 

23.9 

YES 

3.8 X 10-4 

4.3 X 10-4 

0.32 

0.31 

YES 

1.1 X 10-4 

0.88 

A3.0 

35.2 

YES 

2.8 X 10-4 

4 X 10-4 

0.28 

0.29 

YES 

1.3 X 10-3 

0.89 

CA0.7 

27.3 

NO 

- 

- 

- 

- 

YES 

5.9 X 10-4 

0.16 

CA1.5 

29.8 

NO 

- 

- 

- 

- 

YES 

2.6 X 10-3 

0.81 

CA3.0 

17.9 

NO 

- 

- 

- 

- 

YES 

5.5 X 10-3 

0.83 

PE0.7 

24.8 

NO 

- 

- 

- 

- 

NO 

1.5 X 10-4 

0.89 

PE1.5 

25.9 

NO 

- 

- 

- 

- 

NO 

1.9 X 10-3 

0.91 

PE3.0 

31.2 

NO 

- 

- 

- 

- 

YES 

6.7 X 10-3 

0.75 

PF0.7 

27.8 

NO 

- 

- 

- 

- 

NO 

8.6 X 10-3 

0.91 

PF1.5 

27.8 

NO 

- 

- 

- 

- 

NO 

8.8 X 10-4 

0.97 

PF3.0 

21.1 

NO 

- 

- 

- 

- 

YES 

5.6 X 10-3 

0.91 


easy to identify. All their quantities presented have been 
computed using the median values over the last binary 
orbit. We observe that the disc masses are very similar for 
different impact parameters, which means that this quantity 
is -at least in this respect- independent of the amount of 
non-accreted gas. The inner edge of all these mini-discs 
extend from the accretion radius of our SMBH (see Section 
2.2) to around the radius of the Hill sphere^, as expected. 

The circumbinary discs are more difficult to analyse. 
It is not possible to determine at the time we stop the 
simulations which fraction of the gas is going to form a stable 
disc. Hence, we first simply visually establish whether there 
is gas orbiting the binary in closed orbits - this criterion 
is stated in Table 1. Additionally, to estimate the possible 
disc properties, we identify the gas particles that could 
become part of a circumbinary disc. These are going to 
be the particles that are bound to the binary, and have 
orbits calculated around the CoM with pericentre distances 
b larger than a threshold value dehned to be rmin = 2 a for 
the aligned cases and rmin = a for the rest. Gas particles 
with a pericentre distance smaller than this radius will 
either be re-ejected in a slingshot process that prevents 
the gas to complete an orbit around the binary (see e.g. 
top row, left panel of Fig. 3), or accreted, or will become 
part of one of the minidiscs, if present. For ah gas particles 
fulhlling the b > rmin criterion, we measure as before the 
total mass and median eccentricity from the last snapshot 
in each simulation. The amount of gas available to form 
a circumbinary disc increases with the impact parameter, 
simply because gas with larger angular momentum avoids 
being accreted. The median eccentricity of the gas is high 
for most cases due to the initial orbit of the cloud: after the 
hrst passage most of the gas forms a very eccentric tail. The 
exception is case CA0.7, here, due to the small pericentre 
distance, most of the gas is accreted, leading to a very light 
tail respect to the circumbinary ring. The other cases have a 


^ th ~ 0.3a for a circular, equally-massive binary 


more massive tail, which make the eccentricity distribution 
skew towards high values. 

In conclusion, the interaction of individual clouds with 
a binary results in circumbinary discs that are initially very 
eccentric (e > 0 . 8 , note that this estimate includes bound 
material within the highly eccentric tails, the long-term 
evolution of eccentricity will be discussed in Section 6 ) . 
Moreover, by our choice of parameters, the disc masses are 
at most roughly half of the initial cloud mass, and in most 
cases much less than that, only a few percent. This does 
not only make them hard to detect directly, but also implies 
a small effect on the secular evolution of the binary. This 
conclusion is likely to change, however, if we consider the 
cumulative effect of many cloud infah events, or the effect 
of individual massive clouds. The results of that study will 
be presented in a follow-up paper. 


4 ACCRETION RATE AND TOTAL 
ACCRETED MASS 

In Fig. 4 we show the accretion rate and cumulative accreted 
mass for every inclination and impact parameter. Note that 
throughout this section the results are presented in code 
units, thus we can re-scale the results to a range of binary 
masses and periods (for physical rescaling of our results see 
Section 7). 

It is important to emphasise that the accretion 
radius set for each SMBH is very large compared to the 
Schwarzschild radius, so the accretion rate that we compute 
does not correspond to actual accretion onto the SMBHs. 
The gas within our accretion radius will have non-zero 
angular momentum, thus should settle on an accretion disc 
and evolve on a viscous timescale (tvisc), which is typically 
longer than its dynamical timescale (tdyn) in the SMBH 
gravitational potencial. In consequence, our accretion rates 
correspond to the rate at which the gas is added to the 
BH-accretion disc systems, which are unresolved in our 
simulations. For an a-disc (Shakura & Sunyaev 1973), these 


MNRAS 000, 1-17 (2015) 









Formation of discs from infalling clouds onto SMBH binaries 7 


S 0,- 




0.8 

0.5 





-- 0.4 




0.6-^ 

cC 




s 

J 

^ 0.3 




0.4 <1 

o 


m 



0.2 








AO.7 model 

0.2 


ll/Jf 


— black hole #1 

— black hole #2 

— binary 


0.1 


10 12 14 


time [Pq] 



CA0.7 model 

— black hole #1 

— black hole #2 
binary 


10 12 14 16' 


time [Pq] 












. . 








PE0.7 model 




— black hole #1 


I'jf 


— black hole #2 




— binary 


10 12 


time [Pq] 




0.6^ of 0.3 

0.4 < S 0.2 


0.2 0.1 


A1.5 model 

■ black hole #1 

■ black hole #2 

■ binary 


10 12 14 16' 


time [Pq] 




CA1.5 model 




1 


— black hole #1 

— black hole #2 

— binary 

0.8 

0.7 











0.6^ 

iy 

0.4 


1 ^ 



0.4 < 

5 0-3 







P 


0.2 

0.0 

0.2 

0.1 

0.0 




time [Pq 



Figure 4. Evolution of the accretion rate (solid, left axes) and cumulative accreted mass (dashed, right axes) onto the binary (black) 
and each black hole (blue and red), for all simulations. From left to right: pericentre distances of 0.7, 1.5 and 3i7bin- From top to bottom: 
aligned, counter-aligned, perpendicular edge-on and perpendicular face-on orbits. Notice that the range of the left y-axes is not the same 
for all plots. On each curve we see a periodic behaviour that is directly related with the binary period, and in most cases the accretion 
on to both SMBHs are in counter-phase with each other. 


timescales are related as following: 


tvisc 


1 

a 



(4) 


where a is a dimensionless number quantifying the strength 
of the viscosity (with typical values of~ 0.01 — 0.1) and H/R 
is the disc aspect ratio. On the other hand, the dynamical 


timescale (at the sink radius) can be related with the orbital 
period of the binary as following: 



(5) 


where rsink = 0.1a is the accretion radius, M is the binary 
mass and Mi = 0.5M is the mass of one SMBH. 
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8 Goicovic et al. 


As a rough estimate, we would expect the actual 
accretion rate onto each SMBH be related with the binary 
orbital period if the latter is longer than the viscous 
timescale. In other words, when 

«(f) >5xl0-^ (6) 

which implies that there are combinations of parameters 
(e.g., a ~ 0.1 and H/R ^ 0.2) for the unresolved accretion 
discs where the variability presented below would indeed 
affect directly the actual accretion on to the SMBHs (Sesana 
et al. 2012). For systems where the viscous time is longer 
than the binary orbital period, the streaming periodicity 
will not be directly reflected into a variable accretion onto 
each SMBH, and a pair of persistent (unresolved) mini-discs 
will form. However, even in these cases we could expect to 
find observational signatures of this variability (§7). 

Another caveat of our model is that we do not 
include any type of feedback from the SMBHs; radiation 
pressure could reduce the amount of gas that is accreted. 
Nevertheless, the amount of gas that actually reaches the 
SMBHs will likely be a monotonic function of the value that 
crosses our sink radius. 

In summary, the behaviour of the accretion rate shown 
in this section can still be useful to characterise the accretion 
onto the SMBHs and the mini-discs emission. 

As expected, the accretion rate is very high at the 
beginning of the interaction with the cloud. Most of the 
cloud is engulfed by the binary during the first few orbits 
4 — 7), as we can observe in the cumulative mass of each 
figure. This stage corresponds to the first passage of the 
cloud. After that interaction the accretion drops sharply, by 
2 or 3 orders of magnitude in some cases. 

For most orbital configurations there is a clear 
periodicity on the accretion rates. Each black hole has 
an accretion periodicity which matches the binary orbital 
period, but in phase opposition with each other, so that the 
total binary accretion rate has a period of half an orbit. This 
feature is associated with the stream of gas that remains 
bound to the binary and feeds each black hole alternatively 
when they cross the stream. Variability related to the binary 
orbital period is the type of behaviour that we expect will 
help to identify and characterise these systems. 

The accretion in the cases of perpendicular orbits (edge 
and face-on, two lower rows of Fig. 4) tends to be more 
extended than the parallel ones (two upper rows), in the 
sense that there are still significant peaks after the first 
passage of the cloud. This is because, as explained in the 
previous section, the slingshot is less efficient when the 
encounters have perpendicular relative velocities, allowing 
more material to remain around the binary in close orbits. 

For all the orbital configurations we model, we obtain 
two trends with increasing impact parameter: (i) the 
accretion rates and the total accreted mass decrease, and (ii) 
the relative amplitude of the accretion rate peaks during the 
cloud first passage decreases compared to the later stages. 

The compilation of the total mass accreted by the 
binary at the end of each simulation is shown in Fig. 5. Here 
we can see how the accretion is dramatically reduced when 
we increase the impact parameter of the cloud. For instance, 
in the PE3.0 model the total accretion is around 5% of 
the cloud mass, in contrast with the ~ 70% on the smaller 



Figure 5. Fraction of the total cloud mass accreted by the binary 
as a function of the pericentre distance, for the four different 
orbital orientations. Notice that the total mass is strongly 
dependent on the impact parameter of the binary, dropping by 
one order of magnitud in some cases. 

pericentre distance (PE0.7). This is very interesting as it 
shows the transition between a “prompt accretion" regime, 
and one characterised by the formation of circumbinary discs 
for the aligned and perpendicular edge-on orbits (see Section 
6 ). 


5 MISALIGNED MINI-DISCS 

As we mention in Section 3, the only configurations that 
show the formation of extended and stable mini-discs, 
limited by the Hills radius, are the aligned models, which 
is the first case we study below. However, we also measure 
the direction of the gas around the SMBHs for the other 
inclinations to see if possible unresolved mini-discs have 
some preferential orientation. 

5.1 Mini-discs for the aligned configurations 

In all our models with aligned orbits we see the formation 
of prominent and persistent mini-discs around each SMBH 
(top row of Fig. 3). In order to measure the level of 
alignment with the orbit of the binary we compute the 
mini-disc direction using the total angular momentum vector 
of the gas particles within the Hill sphere. We represent 
this direction with the spherical coordinates (^, 0) in the 
reference frame of the corresponding SMBH. The time 
evolution of the mini-disc directions is shown in the Hammer 
projections of Fig. 6. 

From these projections we observe that the mini-discs 
are well defined and evolve smoothly with time, and also 
that they are roughly aligned with the binary orbit (and the 
original cloud orbit), although there is always some level 
of misalignment. Studying closely the time evolution of the 
minidisc orientations, we notice the following behaviour: not 
only the mini-discs precess around the aligned position, but 
on top of that we observe another low-amplitude, periodic 
motion that hereafter we refer to as “wobbling". In order to 
study these motions in more detail, we describe the direction 
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Figure 6. Time evolution of the directions of both mini-discs shown in Hammer projections for the aligned orbits and the three different 
impact parameters. Upper left: O.Ti^binj upper right: 1.5-Rbin cind bottom: 3-Rbin- The circles and squares correspond to the times when 
each mini-disc appears and the end of the simulation, respectively; while the black diamonds correspond to the binary orbit orientation. 
The inset panels show zoom-ins of the projections, in which is noticeable that the evolution of the mini-discs show two combined effects: 
a steady precession around the aligned position and a super-imposed wobbling. 


of the mini-discs with other two angles: inclination and 
position. The first one is the angular difference with respect 
to the aligned position, while the second is the angle between 
the projection of the mini-disc on the orbital plane and an 
arbitrary vector on the same plane. 

Since a similar behaviour is observed in all three runs 
with aligned orbits (see Fig. 6), we now concentrate our 
analysis in the A3.0 model, which ran the longest and has 
the better defined direction evolution. Notice that this run 
also shows a slow but clear alignment of the mini-discs with 
the binary orientation. However, we do not discuss this in 
detail as it depends on the numerics (see below). We show 
the time evolution of the inclination and position angles in 
Fig. 7. We observe that the precession is steady (constant 
slope in absolute value), with a period of around 20 binary 
orbits. On the other hand, the mini-discs inclination tend 
to decrease with time and to wobble with a period almost 
exactly half of the binary period. As we show below, this 
behaviour is expected for misaligned discs. 

5.2 Dynamics of misaligned discs 

Based on the study by Bate et al. (2000) we can obtain an 
analytical understanding of the dynamics of our misaligned 
mini-discs around each SMBH. We consider a circular binary 
with components Mi and M 2 , with separation of a. We sit 
on a non-rotating frame with z axis parallel to the rotation 
axis of the binary, centred on one of the SMBHs. The force 
at a position vector f (with r a) caused by the other black 



time [Pbin] 


Figure 7. Upper panel: inclination angle evolution of the 
mini-discs respect to the aligned position. Both the slow decline 
due to dissipation and the wobbling discussed in the main text 
are clearly visible. Lower panel: position angle evolution, i.e. 
angle between the projection of the mini-disc direction on the 
orbital plane and an arbitrary vector on the same plane. This 
evolution shows how the discs slowly and steadily precess around 
the aligned position. The period of this movement is around 20 
binary periods. 

hole is given by: 

- GM2mr ^ ^ 3GM2mr.^ 

F 2 =-^- r-\ -r- Ir • a)a, (7) 
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where a is the position vector of the secondary SMBH and 
rrir is a test mass. 

The right side of equation (7) can be split into 
2 contributions: the isotropic term {m — 0) and the 
quadrupole (m — 2), and each of them can be associated 
with a different effect. 


(i) The m = 0: If we consider a ring of gas with mass rrir 
and radius Ur, the torque produced by this term is: 


^“ = -4V 


3 ( GM2mral 


sin (5 cos 5i. 


(8) 


The angular frequency of a ring of disc material is given 


by 


Q^d — —Q^d sin 5j + Q^d cos 8k^ 


(9) 


where Q^d is a function of the radius of the considered ring, 
5 is the inclination angle of the ring respect to the binary 
plane. Then, as - To = 0, the net effect of this torque is to 
produce a precession around the z axis. The mean precession 
rate is 

3 ^ /ar\^ 

_ = _,cos5(-) , (10) 

where q is the mass ratio of the binary. 

(ii) The m — 2\ In this case the torque on the ring is 

tL = — ^ A ^ ^ cos{2Qbt)t 

4\ J (11) 

+ cos S sm{2Qbt)j + sin S sin(2Qbt)^], 

where is the binary angular frequency. This oscillating 
torque is also perpendicular to the rotation of the ring (T 2 • 
= 0), and the effect is to produce an oscillation around 
the steady precession, with a frequency equal to twice the 
binary frequency, consistent with what we measure for the 
inclination in Fig. 7. The amplitude of the wobble is roughly 
~ Up/ (2Qb)- 


The behaviour of the entire disc will be given by the 
integral of each ring. Then, the net precession rate is 


g=i^cos5g(^) , 

where 

3 0/0 

K= - 

4 ff Sr-3/^dr 


( 12 ) 


(13) 


and E is the mass surface density of the disc. Some typical 
values of K are: iF = 15/32 0.487 for a constant surface 

density profile (Larwood et al. 1996); and K = 3/10 for 
E (X (Hartmann et al. 1998). 

Thus typically 


Up 

Qb 




3/2 

cos (5 


(14) 


Evaluating equation (14) with the approximate values 
of our mini-discs we find 



and also the amplitude of the wobble in the outer part of 


the disc, which is the dominant contribution, is 

^wob ~ « 0.02 rad « 1 deg. (16) 

Both quantities are actually very close to what we observe 
on Fig. 7. 

All these calculations are made assuming that the discs 
are able to communicate the precession efficiently without 
breaking, which is clearly the case of our simulations as we 
observe them moving as whole. However, the ability of the 
disc to process rigidly depends strongly on its aspect ratio 
and viscosity (Papaloizou & Terquem 1995; Larwood et al. 
1996; Lubow & Ogilvie 2000; Fragner & Nelson 2010; Dogan 
et al. 2015), which are not well resolved quantities in our 
models due to the small number of particles (~ 1000) that 
shape our mini-discs and the numerics itself. In addition, the 
presence of the accretion radius around the sink particles 
excises the inner portion of the mini-discs, which does not 
allow US to model the precession and wobbling particularly 
onto those scales. 

The numerical models performed by Fragner & Nelson 
(2010), study the evolution of individual discs that arise 
in misaligned binary systems, showing precession and 
wobbling (i.e. periodic perturbation in the inclination 
angle), similarly as we found in our mini-discs. Their 
grid-based hydrodynamic code allows them to control 
better the viscosity’s influence in the disc evolution. 
They show that, for several combinations of aspect ratios 
and viscosities, the discs will efficiently communicate the 
differential precession and move as a rigid body. In 
particular, they find that thin discs (h < 0.03) with high 
viscosities will achieve a state of rigid precession after 
developing a twist. More important, when the discs are not 
disrupted due to the differential precession, the periodic 
perturbations in the inclination are able to travel all the 
way to the central parts (cf. their Fig.6), which we are not 
able to resolve. 

In conclusion, even if we are not able to robustly study 
the global dynamics of the mini-discs due to the relatively 
low resolution we can afford with our models, we expect that 
the misalignment arising from the infall of extended portions 
of gas will produce the precession and wobbling we observe. 


5.3 Mini-discs for other inclinations 

As explained in Section 3, during the whole duration of 
the simulations with different inclinations (counter-aligned 
and perpendiculars) we do not observe the formation of 
persistent mini-discs, probably because we do not have 
the spatial resolution due to the artificially large accretion 
radii and the low angular momentum of the captured gas 
relative to the individual SMBHs. Nevertheless, we can still 
measure whether the material around each SMBH has some 
preferential direction. We do so by computing the angular 
momentum of all accreted gas between outputs, which is 
an indicator of the direction that the unresolved mini-discs 
might have, if they exist. We consider the time after the first 
passage of the cloud (~ 8 orbits) for every snapshot until 
the simulation ends. We show an example of the projections 
obtained with the larger impact parameter in Fig. 8, where 
each point represents the direction of the gas in a particular 
output. 
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Figure 8. Evolution of the direction of all the accreted gas by each SMBH (red and blue dots) in a Hammer projection. The different 
plots correspond to different inclinations: counter-aligned (upper left), perpendicular edge-on (upper right) and perpendicular face-on 
(bottom). As indicated in the legend, all correspond to the larger impact parameter. The black diamonds indicate the binary direction. 
Notice how in all cases the gas tends to follow the initial orientation of the cloud, located always at (90°, 0°), and not the orientation of 
the binary. The effect is weakest for the counter-aligned orbit case. 


For the counter-aligned orbits (upper left panel of 
Fig. 8) the different orientations of the gas are not 
concentrated around any particular point. However, they 
are preferentially on the right side of the projection, 
which indicates that the gas tends to be counter-aligned 
respect to the binary rotation. On the other hand, for 
both perpendicular orbits (upper right and bottom panel of 
Fig. 8) the gas tends to cluster more clearly around (90°, 0°), 
which is the initial direction of the cloud orbit. This is 
because the dynamical interaction with the binary is not 
able to efficiently change the angular momentum direction of 
the surrounding gas on such short time-scales. In conclusion, 
the mini-discs that might arise from these perpendicular 
accretion events will be completely misaligned respect to 
the binary orbit, but they are likely to be roughly aligned 
with each other. 

Due to the large misalignment between the possible 
mini-discs and the binary, we expect other periodic 
effects to appear, like the Kozai-Lidov oscillations (Kozai 
1962; Lidov 1962) for example, where a test particle 
around one component of the binary periodically exchanges 
its inclination for eccentricity. Using hydrodynamical 
simulations, Martin et al. (2014) showed that this effect 
is also present in fluid discs. The time-scale for these 
oscillations is expected to be several times the binary orbital 
period - Martin et al. estimate txL ~ 17(0.35a/i?out)^'^^U6 
for an equally-massive, circular binary and a mini-disc with 
surface density given by S oc which implies periods 

over 20 binary periods for mini-discs inside the Hill radius. 


However, the changes on the disc inclination and eccentricity 
are large, which could have implications on processes such 
as the shaping of jets, the feeding onto the SMBHs and star 
formation. 

In summary, as hinted by our set of simulations, the 
misalignment appears to be a natural outcome from infalling 
cloud events, even when their orbits are aligned with 
the binary. Through the particular dynamics due to the 
interaction with the non-Keplerian gravitational potential 
of the binary, this could have important implications on the 
observability of these systems (see Section 7). 


6 CIRCUMBINARY DISCS 

As we showed in Section 3, a circumbinary disk promptly 
forms ^ whenever the impact parameter of the infalling 
cloud is large enough (or when its orbit is retrograde 
to that of the binary). However, our simulations are too 
short to assess the physical properties of these discs as 
they evolve toward a (possibly) steady state. In order to 
investigate the longer term evolution of these discs, we take 
two representative cases (namely, cases A3.0 and PE3.0 in 
table 1). We re-simulate them for 30 further BHB orbits, 
using the final snapshots as initial conditions, but keeping 
only a sub-set of particles, to save computational time. We 
selected only particles bound to the binary with a period 

^ within the 20 orbits or so that we are modelling 
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Figure 9. Column density maps of the PE3.0 re-simulation. 
The different times correspond to the beginning and the end of 
the re-simulation. The upper and lower panels show face-on and 
edge-on views of the disc, respectively. 


smaller than 30 times the binary period, as the contribution 
of the excised particles is dynamically negligible. 

The PE3.0 model is particularly interesting, because the 
circumbinary disc retains memory of the original orientation 
of the cloud and it is, therefore, perpendicular to the binary 
orbit. Column density maps at the beginning and the end of 
the re-simulation, from two edge-on views of the binary, are 
shown in Fig. 9. In the lower panel, where the disc is seen 
roughly edge-on, we observe how it tends to slowly align 
to the binary orbit. In order to measure this evolution, we 
use the total angular momentum of the gas within a fixed 



Figure 10. Time evolution of the gas inclination on the 
re-simulation of the PE3.0 model. Because the circumbinary disc 
is not a well-defined structure, we sum the angular momentum 
of all the gas, and also only that within two different radii, as 
indicated on the legend. 


radius as a proxy for the direction of the disc. The evolution 
of the inclination angle is shown in Fig. 10 for different 
definitions of the disc extent, showing that this result is 
robust. We estimate that the alignment timescale is around 
1000 orbits, should the evolution remain roughly constant. 
However, as explained for the mini-discs evolution, this 
timescale will depend critically on the viscosity treatment, 
and its study is beyond the scope of this paper. We 
also notice that the inclination evolution shows the same 
oscillations that the mini-discs showed in the aligned case. 
These oscillations have also half of the binary period, as 
expected since the dynamics is driven by similar processes 
as in the former case. This periodic perturbation has been 
overlooked on most studies of misaligned discs because it 
does not have a secular effect on its evolution. However, 
our model suggests that it might have some interesting 
implications. For example, the density waves produced by 
the oscillatory perturbations enhance fragmentation on the 
gas, that we can see in the form of clumps in Fig. 10. Some 
of these clumps may form stars, a fraction of which will 
end up producing observable tidal disruption events (TDFs; 
see Section 7). Another possible signature of this oscillation 
might be imprinted in the shifting of spectral lines, specially 
coming from the inner regions where their amplitude is 
larger. The longer simulation also allow to investigate trends 
in the circumbinary disc eccentricity. By inspecting the 
face-on views (upper panels of Fig. 9) it seems that the 
material is more concentrated on eccentric orbits towards 
the end of the simulation. In order to measure this, we plot 
the eccentricity distribution at three selected times in the 
upper panel of Fig. II. The distribution becomes narrower 
as the simulation advances, but actually keeps essentially 
the same median eccentricity of ~ 0.6. 

For the re-simulated A3.0 model, column density maps 
for the beginning and the end are shown in Fig. 12. As 
expected, the gas and the binary are essentially coplanar, 
with an initial inclination of only 4°, which decreases by 
~ 0.5° during the re-simulation. The eccentricity evolution, 
shown in the lower panel of Fig. II, is more interesting. Here 
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Figure 11. Eccentricity distribution of the gas for PE3.0 
(top) and A3.0 (bottom) re-simulations at different times. Each 
distribution is normalised by its maximum value. The vertical, 
dashed lines indicate the median of each distribution. For the 
perpendicular disc we observe that the distribution narrows, 
keeping roughly constant its median value, while for the aligned 
disc the distribution shifts towards higher values with time. 



Figure 13. Eccentricity distribution of the gas for the 
counter-aligned simulations at different times, increasing impact 
parameter from top to bottom. Each distribution is normalised by 
its maximum value. The vertical, dashed lines indicate the median 
of each distribution. In the three models we observe the material 
decreasing its eccentricity with time, reaching values close to zero. 
This process occurs faster for smaller impact parameters. 



of Fig. 3). In this configuration, the formation of a 
clear circumbinary structure is extremely quick, and mass 
adds-up as the very eccentric tail joins it. We show the 
evolution of the gas eccentricity distributions in Fig. 13. At 
the beginning the distribution is completely skewed towards 
high eccentricities due to the initial conditions, but after 
the interaction with the binary the gas eccentricity shifts 
towards lower values, and a striking bi-modality appears. 
This is clearly related to the formation of the ring, and it 
occurs on shorter time-scales for smaller impact parameters. 

All our circumbinary discs appear to evolve differently 
according to their relative inclination to the binary. In the 
aligned model, the gas tends to increase its eccentricity; 
in the counter-aligned it tends to become circular, while 
in the perpendicular the eccentricity retains its value. 
This is driven by the dynamical interaction with the 
binary: a prograde encounter of a gas particle with one 
of the SMBHs will tend to increase its specific energy and 
angular momentum, increasing the eccentricity; on the other 
hand, a retrograde encounter will work in the opposite 
direction, circularising the material. For the perpendicular 
encounters, the binary is unable to change the orbital 
angular momentum of the gas, keeping its eccentricity 
roughly constant during the interaction. 


Figure 12. Column density map at the beginning and the end of 
the A3.0 re-simulation. In this case we only show a face-on view 
of the disc and the binary, because they are almost completely 
aligned. Here we see that the disc increases its eccentricity. 

we observe a different behaviour respect to the model PE3.0; 
the orbits become more eccentric as time advances, which 
is reflected in the shifting of the distribution towards higher 
values and a clear increase of the median. 

Finally, a clear circumbinary ring also appears in 
our simulations of the counter-aligned cases, for all the 
investigated initial parameters (see upper middle row 


7 PHYSICAL SCALING AND 

OBSERVATIONAL CONSEQUENCES 

The efforts to confirm observationally the existence of 
SMBH binaries have increased in the last few years, 
motivated by their key role as probes of the hierarchical 
growth of galaxies (Sesana et al. 2011). Observing 
these objects is unfortunately very challenging, as their 
separations cannot be resolved by current capabilities on 
most galactic nuclei. Additionally, for the very few existent 
candidates, the observed signatures can also be explained 
by alternative scenarios that do not require a BHB (for a 
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compilation of candidates and prospects in observational 
searches for BHBs, see Dotti et al. 2012 and Bogdanovic 
2015). 

In the late stages of their evolution, when the binary 
separation is sufficiently low (< 10“^ pc), gravitational 
waves will efficiently extract its angular momentum and 
energy, rapidly leading it to coalescence. The enormous 
amount of gravitational wave emission will be detectable 
by pulsar timing arrays (Sesana 2013) or by space-based 
missions like eLISA (Amaro-Seoane et al. 2013a). However, 
these detections will need to be complemented with 
observations in the electromagnetic spectrum in order to 
localise and characterise the sources. 


7.1 Physical scaling 

In order to link our results with possible observational 
signatures we scale all physical quantities from code to 
physical units. To perform the scaling, we have to choose two 
parameters: the mass of the binary and the critical density 
for the equation of state (eq. 3); fixing these two values will 
determine the rest of the units. 

In Table 2 we present physical scaling down to 84 years. 
Shorter periods would require too large densities compared 
to what is observed for molecular clouds. Moreover, for 
the more massive systems, a period of 84 years already 
implies a cloud temperature of ~ 10^ K, which is already 
unrealistically high (see e.g. Meijerink et al. 2007) and 
cannot be pushed further. A period of roughly a century 
is too long to look for the variability associated to it 
in observed data. Nevertheless, even if our model is 
unable to represent directly more compact systems, we 
still expect that the behaviour we find is qualitatively 
representative of those more rapidly varying systems, and 
below we discuss several possible observational signatures. 
In upcoming publications we will present models designed 
specifically for more compact systems. 


7.2 Observational signatures 

The feeding rates onto the binary and each SMBH show 
variability for all the configurations we model, always related 
with the orbital period. Recall that we measure the accretion 
rates at the sink radius Rsink, which is large compared 
to the Schwarzschild radius. Still, what AGN observations 
reveal is the luminosity of the accretion disc at different 
radii depending on the measured wavelength. We expect the 
variable accretion rate we obtain at Rsink to represent actual 
variability of the accretion disc at large radii, which in some 
AGN shows up in optical or infrared light curves (e.g.. Lira 
et al. 2011). 

Variability due to binary feeding could thus appear in 
AGN light curves (see Sillanpaa et al. 1988, for the iconic 
case of OJ 287), and be detected with future time-domain 
surveys like the Large Synoptic Survey Telescope (LSST; 
Ivezic et al. 2008). Interestingly, periodicity has been 
recently claimed for a couple of systems proposed as binary 
candidates (Graham et al. 2015; Liu et al. 2015). 

The presence of the mini-discs in our simulations, with 
their misalignment and evolution, allows us to explore novel 
and promising observational features. Graham et al. (2015) 


showed that the blazar PG 1302-102 has a roughly sinusoidal 
light curve with a period of 5 years, and mentioned 
as a possible explanation the precession of a jet. If we 
assume that the orientation of a possible jet is given by the 
mini-disc®, we expect to observe variability related to the 
mini-disc wobbling found in our models and predicted earlier 
by Bate et al. (2000), specially if the jet is close to the line 
of sight. Graham et al. (2015) fitted the observed light curve 
with a wobbling amplitude of around 0.5° which is actually 
very close to what we find in Fig. 7. Therefore, we suggest 
that the observed variability of this source could be due to 
mini-disc wobbling as shown in our models. Notice that this 
interpretation would imply a rest-frame orbital period for 
the binary of eight years, one order of magnitude longer than 
under the assumption that the variability is due to lumps in 
the circumbinary disc (D’Orazio et al. 2015). 

In the cases where the viscosity of the individual discs 
is low, the differential precession induced by the companion 
SMBH is not communicated efficiently throughout the discs, 
which then could break. The disruption of the disc into 
rings that process independently could also be a source 
of variability, as the dissipation between the gaseous rings 
is enhanced and this promotes stronger accretion onto 
the central SMBH (Dogan et al. 2015). There is recent 
observational evidence for this process in a proto-planetary 
disc (Gasassus et al. 2015). 

Additional observable features could be due to warps 
in the misaligned mini-discs. We do not find warps in our 
simulations due to the relatively low resolution we can 
afford. However, they are expected due to the gravitational 
pull of the companion (Moeckel V Bally 2006). The presence 
of a warp can cause part of the disc to block other parts of 
the disc or a central source. Additionally, it can change the 
viewing angle which is important for synchrotron emission. 
Both effects occur periodically, on the dynamical time-scale 
of the disc. The effects of warps have been observed in 
several astrophysical objects such as active galactic nuclei, 
circumstellar discs and stellar mass black holes through 
variability in (i) the photometry (e.g. Herrnstein et al. 2005; 
Manset et al. 2009; Bouvier et al. 2013); (ii) the spectrum 
(e.g. Reynolds et al. 2009; Looper et al. 2010); and (hi) 
polarisation (e.g. Manset et al. 2009; Roland et al. 2009; 
Gheng et al. 2015). 

Gircumbinary discs can also have variable emission. In 
particular. Figs. 9 and 12 show a series of discrete shock 
fronts propagating from the binary into the tail of the 
eccentric disc. Such regular feature is not seen in comparable 
simulations featuring a single SMBH (Bonnell & Rice 
2008). Material approaching the periastron is accelerated 
and flung away into the tail when it is in phase with 
one of the two SMBHs. This creates ejected waves with a 
periodicity of half the binary orbital period that compress 
and shock the surrounding material, possibly leading to 
periodic enhancements in luminosity and discrete episodes 
of star formation. As the gas is located farther from the 
black holes (r > 2a, for a prograde disc), we expect the 
period of such variability to be typically longer than the 


® The jet orientation could also be given by the spin of the SMBH, 
which is not necessarily related with the mini-disc orientation, 
specially at these scales. 
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Table 2. Compilation of physical units of the initial values on our simulations. From left to right: binary mass, critical density for EoS, 
separation of the binary orbit, separation in terms of Schwarzschild radius, binary period, cloud mass, initial distance between the cloud 
and the binary, cloud radius, modulus of the cloud initial velocity, cloud velocity at the periastron for the smaller impact parameter, 
initial temperature of the cloud. 


Mbin (Mq) 

Per (gem ^) 

a (pc) 

c /Rsch 

Bbin (yrs) 

Mci (Mo) 

Rcl (pc) 

'Cini (km/s) 

'Cperi,0.7 (km/s) 

Tini (K) 

10® 

10-14 

0.2 

10® 

8370 

O 

1 — 1 

0.5 

40.8 
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100 

10® 
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0.04 

2 X 10® 
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10^ 

0.1 

85.7 

720 

470 

10® 

10-10 

0.009 

5 X 10"! 

84 

O 

1 — 1 

0.02 

187.7 

1577.3 

2170 

1- 

O 

1 — 1 

O 

1 

0.4 

2 X 10® 

8370 

10® 

1 

85.7 

720 

470 

10'^ 

10-12 

0.09 

5 X 10"! 

837 

10® 

0.2 

187.7 

1577.3 

2170 

1- 

o 

1 — 1 

O 

1 

o 

0.02 

10^ 

84 

10® 

0.05 

408 

3429 

10^ 

10® 

10-14 

0.9 

5 X 10"! 

8370 

10® 

2.3 

187.7 

1577.3 

2170 

10® 

10-12 

0.2 

10^ 

837 

10® 

0.5 

408 

3429 

10^ 

10® 

10-10 

0.04 

2 X 10® 

84 

10® 

0.1 

877.2 

7372.5 

4.6 X 10^ 


binary period (e.g. Farris et al. 2015). Then, a source 
with a photometric period shorter than the spectroscopic 
period could be interpreted as a binary surrounded by 
a disc - the binary varies on one, or half, an orbital 
period due to accretion, while the circumbinary varies on 
its own dynamical time. Moreover, the ratio between both 
variability periods could be a tool to study the geometry 
and extension of the gas around the binary. 

We have shown that increasing the impact parameter 
results in more clumping for the gas. The presence of clumps 
could have indirect implications on the observability of a 
binary through star formation and posterior tidal disruption 
events (TDEs). Amaro-Seoane et al. (2013b) showed that 
fragmentation in a circumbinary disc results on in-situ star 
formation and an increase of the rate of TDEs respect 
to what is expected for typical galaxies. Later on, Brem 
et al. (2014) demonstrated that the presence of the binary 
instead of a single SMBH will produce a distinctive signature 
on the reverberation mapping after such an event. Based 
on this, we might expect to detect these circumbinary 
structures through TDEs. Although the rate will depend on 
the efficiency of star formation on the gas tail and the disc 
itself, which we do not model here. In any case, LSST will 
detect thousands of TDEs, greatly increasing the chances of 
detecting this kind of events even if they are relatively rare. 

Finally, star formation in a circumbinary disc could 
also refill the loss-cone and affect the evolution of the 
binary orbit by exchange of energy and angular momentum 
with the stars via 3-body interactions (Amaro-Seoane et al. 
2013b). Moreover, eccentric discs of stars, which is the most 
likely output of the near-radial infall of clouds that we are 
modelling (Figs. 11 and 13), will be subject to instabilities 
(see e.g. Madigan et al. 2009) that can increase even further 
the amount of stars with orbits that will interact directly 
with the binary, enhancing the probability of TDEs and/or 
evolving the binary orbit. 


8 SUMMARY 

We presented numerical, SPH models of the evolution 
of turbulent clouds in near-radial infall onto equal-mass, 
circular supermassive black hole binaries. In order to explore 
different orbital configurations for the cloud we performed a 


total of 12 simulations, changing both the impact parameter 
and the relative inclination between orbits. 

We studied the formation of discs around the binary 
and each SMBH depending on those orbital configurations. 
Our main findings can be summarised as follows: 

(i) We only observe the formation of stable and prominent 
(up to the Hill radius) mini-discs for the aligned models, 
independent of the cloud impact parameter. For the other 
orientations, where the gas trapped by the individual 
SMBHs has little angular momentum, we do not have 
the spatial resolution to observe the possible formation of 
smaller mini-discs. 

(ii) The misalignment of the mini-discs around each 
SMBH seems to be a natural outcome of the infall of 
extended gas clouds. We have confirmed the analytical result 
that misaligned mini-discs will precess and wobble due to the 
presence of a companion SMBH. 

(iii) When the impact parameter is large enough 
(pericentre distance 3Rbin) we always observe the 
formation of a circumbinary disc, independent of the orbit 
orientation. In the counter-aligned cases there is formation 
of a circumbinary disc (or ring) for all impact parameters. 

(iv) The circumbinary discs tend to follow the initial 
orientation of the cloud. We have confirmed that a 
misaligned disc will evolve towards alignment. 

(v) The circumbinary discs are initially very eccentric 
(e > 0.6), and their early evolution will depend on the 
relative orbital inclination. In the aligned model, the gas 
tends to increase its eccentricity, while the opposite is true 
in the counter-aligned model. In the perpendicular cases the 
eccentricity distribution becomes narrower without changing 
its median value. 

(vi) The feeding rates onto the binary and each SMBH 
show variability for all our models, always related to the 
orbital period. 

Although our simulations do not allow a direct scaling 
to BHBs with orbital periods of few years (i.e., accessible to 
future time-domain surveys), we argue that the observed 
qualitative phenomenology, which is mostly driven by 
the gravitational torques exerted by the BHB onto the 
gas, can be safely extrapolated to shorter periods. Such 
phenomenology has distinctive electromagnetic signatures, 
which are very relevant for the future identification and 
characterisation of BHBs. 
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